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Abstract 

A fast and exact algorithm is developed for the spin ±2 spherical harmonics trans- 
forms on equi-angular pixelizations on the sphere. It is based on the Driscoll and 
Healy fast scalar spherical harmonics transform. The theoretical exactness of the 
transform relies on a sampling theorem. The associated asymptotic complexity is of 
order Oil? log^ L), where 2L stands for the square-root of the number of sampling 
points on the sphere, also setting a band limit L for the spin ±2 functions consid- 
ered. The algorithm is presented as an alternative to existing fast algorithms with 
an asymptotic complexity of order 0(L S ) on other pixelizations. We also illustrate 
these generic developments through their application in cosmology, for the analysis 
of the cosmic microwave background (CMB) polarization data. 

Key words: computational methods, data analysis, cosmology, cosmic microwave 
background 



1 Introduction 

In the last few years, the analysis of the temperature anisotropies of the cos- 
mic microwave background (CMB), together with other cosmological observa- 
tions, has allowed the definition of a precise concordance cosmological model. 
These observations culminated with the release of the three-year data of the 
Wilkinson Microwave Anisotropy Probe (WMAP) satellite experiment. The 

* Corresponding author. 

Email address: yves.wiaux@epfl.ch (Y. Wiaux). 



Preprint submitted to Elsevier 



February 4, 2008 



cosmological parameters are now determined with an unprecedented preci- 
sion of the order of several percent jT82l3P] . In the concordance model, the 
CMB originates from quantum energy fluctuations defined in a primordial era 
of inflation. These tiny fluctuations are Gaussian in first approximation. The 
cosmological principle of homogeneity and isotropy of the universe is also as- 
sumed. The observed radiation is therefore understood as a unique realization 
of a Gaussian and stationary (i.e. homogeneous and isotropic) random pro- 
cess on the sphere, which may be completely characterized from its two-point 
correlation functions, or the corresponding angular power spectra. 

The present concordance values of the cosmological parameters are obtained 
through a best fit of the theoretical temperature angular power spectrum of 
the CMB with the experimental data. Beyond temperature anisotropies, i.e. 
intensity anisotropies, a polarization of the CMB is also present which consti- 
tutes a complementary source of information for cosmology. This polarization 
is produced through Thomson scattering at the epoch of recombination. The 
degree of polarization of the CMB is expected to be of the order of 10 percent 
of the temperature anisotropies at small scales, and lower at large scales. As 
Thomson scattering only produces linearly polarized light, the CMB radiation 
is completely described by its temperature T, and its linear polarization Stokes 
parameters Q and U [5.6.7.8,9. 10J . First polarization measurements were re- 
cently obtained, notably by the WMAP experiment [11] . Future CMB exper- 
iments such as the Planck surveyor satellite experiment will allow a deeper 
probe of the temperature and polarization spectra, thanks to improved sensi- 
tivity and resolution on the whole sky. 

From the mathematical point of view, the observable temperature T is a scalar 
function on the sphere, i.e. invariant under local rotations in the plane tan- 
gent to the sphere at each point. The associated invariant TT angular power 
spectrum results from the decomposition of the temperature in scalar spher- 
ical harmonics. But the observable polarization Stokes parameters Q and U 
transform as the components of a transverse, symmetric, and traceless rank 2 
tensor under local rotations. However, scalar electric E and magnetic B polar- 
ization components may equivalently be defined from the parameters Q and 
U. The associated invariant EE and BB polarization angular power spectra, 
and the cross-correlation TE spectrum result from the decomposition of the 
combinations Q ± iU in spin ±2 spherical harmonics on the sphere (6] . From 
the numerical point of view, the asymptotic complexity associated with a naive 
quadrature based on the definition of the scalar and spin ±2 spherical harmon- 
ics transforms is of order 0(L 4 ), where L roughly identifies the square-root 
of the number of sampling points on the sphere. Corresponding computation 
times for the analysis of megapixels all-sky maps such as those of the ongo- 
ing WMAP or the forthcoming Planck experiments are of the order of days. 
Fast and precise computation methods for the scalar and spin ±2 spherical 
harmonics transforms of functions on the sphere are therefore needed. 
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Beyond cosmology, an algorithm for the spin ±2 spherical harmonics trans- 
forms will find application in the spectral analysis of arbitrary spin ±2 signals 
on the sphere, components of transverse, symmetric, and traceless rank 2 ten- 
sor fields under local rotations. 

In the present work we develop a fast algorithm for the spin ±2 spherical 
harmonics transforms of band-limited functions on the sphere. It is based 
on an existing fast algorithm for the scalar spherical harmonics transform. 
It is defined on 2L x 2L equi-angular pixelizations in spherical coordinates 
(9, ip) on the sphere. The algorithm is theoretically exact thanks to the exis- 
tence of a sampling theorem. The associated asymptotic complexity is of order 
0(L 2 log 2 . L). Corresponding computation times for megapixels maps are re- 
duced to seconds. The algorithm is presented as an alternative to existing fast 
algorithms with an asymptotic complexity of order 0(L 3 ) on other pixeliza- 
tions which are widely used in the context of astrophysics and cosmology. 

In § [21 we recall the notion of spin n functions on the sphere. In § [3], we define 
and implement a fast and exact algorithm with complexity 0{L 2 log 2 L) for 
the spin ±2 spherical harmonics transforms on equi-angular pixelizations. In 
§ HI we illustrate the interest of our algorithm in the context of the analysis 
of CMB polarization data. We finally briefly conclude in § [51 



2 Spin n functions on the sphere 

In this section, we discuss standard harmonic analysis on the sphere and on 
the rotation group SO (3). We also discuss the notion of spin n functions on the 
sphere and their decomposition in a basis of spin-weighted spherical harmonics 
of spin n. 



2. 1 Standard harmonic analysis 

Let the function G(u) be a square-integrable function in L 2 (S 2 , dfl) on the unit 
sphere S 2 . The spherical coordinates of a point on the unit sphere, defined in 
the right-handed Cartesian coordinate system (o, ox, oy, oz) centered on the 
sphere, read as u = (9, if). The angle 9 G [0, tt] is the polar angle, or co-latitude. 
The angle <p G [0, 2ir[ is the azimuthal angle, or longitude. The invariant 
measure on the sphere reads dQ = dcos9dip. The standard scalar spherical 
harmonics Yi m (u), with I G N, m G Z, and \m\ < I, form an orthonormal basis 
for the decomposition of functions in L 2 (S 2 ,dQ) on the sphere [12]. They 
are explicitly given in a factorized form in terms of the associated Legendre 
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polynomials P{ n (cos9) and the complex exponentials e m<p as 



Y lm (9,p) 



21 + 1 (l-m)\ 
47r (l + m)\ 



1/2 



Pr (cos 9) e imv . (i; 



This corresponds to the choice of Condon-Shortley phase (— l) m for the spher- 
ical harmonics, ensuring the relation (— l) m Y^(c<j) = F^_ m )(o;). This phase is 
here included in the definition of the associated Legendre polynomials |T3]fT2] . 
Another convention [T3] explicitly transfers it to the spherical harmonics. Any 
function G(uS) on the sphere is thus uniquely given as a linear combination 
of scalar spherical harmonics: G{uS) = J2i£N J2\ m \<l Gi m Yi m (u) (inverse trans- 
form), for the scalar spherical harmonics coefficients G\ m = J S 2 dfl Y^(u;)(jr(u;) 
(direct transform), with \m\ < I. 

Let now G(p) be a square-integrable function in L 2 (SO(3),dp) on the group 
5*0(3) of three-dimensional rotations. Any rotation p G SO (3) may be ex- 
plicitly given in the Euler angles parametrization as p = (p,9,x), describing 
successive rotations by the Euler angles x G [0, 2tt[, 9 G [0, n], and ip G [0, 2ir[, 
around the axes of coordinate oz, oy, and oz respectively. The invariant mea- 
sure on the rotation group reads dp = dipd cos 9d\. The Wigner D-functions 
D l mn (p), with I G N, m, n G Z, and \m\, \n\ < I, are the matrix elements of the 
irreducible unitary representations of weight I of the rotation group 5*0(3), 
in L 2 (50(3), dp). By the Peter- Weyl theorem on compact groups, the matrix 
elements D l ^ n also form an orthogonal basis in L 2 (SO(3),dp) [12]. They are 
explicitly given in a factorized form in terms of the real Wigner (^-functions 
d l mn (9) and the complex exponentials e~ imip and e~ inx as 

DL(¥>Ax)=e- im *d l mn (9)e-™x (2) 

Any function G(p) in L 2 (50(3), dp) is thus uniquely given as a linear combi- 
nation of Wigner ^-functions : G(p) = E ZeN (2Z + l)/8vr 2 E\ m \,\n\<i ^nC(z') 
(inverse transform), with \m\, \n\ < I and where G l mn = J so ^ dp D l mn (p)G(p) 
(direct transform) stands for the with Wigner D-functions coefficients. 



2.2 Spin n functions 



Let us define a spin n square-integrable function n G(u) in L 2 (S 2 ,dfl) on the 
sphere. The Euler angles (ip, 9, x) associated with a general rotation p in three 
dimensions may also be interpreted in the reverse order as successive rotations 
by <p around oz, 9 around oy', and x around oz", where the axes oy' = oy'((p) 
and oz" = oz"(ip, 9) are respectively obtained by the first and second rotations 
of the coordinate system by ip and 9 [H]. The local rotations of the basis 
vectors in the plane tangent to the sphere at u = (9, ip) are rotations around 
oz", therefore associated with the third Euler angle x- Spin n functions on the 
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sphere n G(u), with n G Z, are defined relatively to their behavior under the 
corresponding right-handed rotations by Xo as [1511161117] : 

n G' (u) = e~ inxo n G (u) . (3) 

The standard square-integrable functions on the sphere considered above are 
spin or scalar functions. Let us emphasize that the rotations considered are 
local transformations on the sphere around the axis oz" = oz"((p, 8), affecting 
the coordinate x i n the tangent plane independently at each point uj = (8, ip), 
and according to x' = X ~ Xo- They are to be clearly distinguished from the 
global rotations by x around oz associated with the alternative Euler angles 
interpretation, which affect the coordinates of the points uo = (8, ip) on the 
sphere. Our sign convention in the exponential is coherent with the definition 
PJ) below for the spin-weighted spherical harmonics of spin n. It is opposite 
to the original definition [15] , while equivalent to recent notations used in the 
context of the CMB analysis [6110] . 

Recalling the factorized form (j2j), spin functions are equivalently defined as the 
evaluation at x = of any function in L 2 (SO(3), dp) resulting from an expan- 
sion for fixed index n in the Wigner D-functions D mn (ip, 8, x)- The functions 
D l mn (ip, 8, 0) or D l £,_ n Jip, 0) thus naturally define for each n an orthogonal 
basis for the expansion of spin n functions in L 2 (S 2 , dfl) on the sphere. After 
normalization in L 2 (S 2 ,dQ), the spin-weighted spherical harmonics of spin n 
are given in a factorized form in terms of the real Wigner (^-functions d l mn (8) 
and the complex exponentials e tmip as 



n Y lm (8, <p) = (-1)" ^ ^d l m{ _ n) (8) , (4) 

with / G N, I > \n\, and m G Z, |m| < I. In particular, the symme- 
try properties of the Wigner ci-functions [12J imply the generalized symme- 
try relation (— \) n+m n Y^(u;) = _„yj(_ m ) (u) . The spin spherical harmon- 
ics explicitly identify with the standard scalar spherical harmonics for the 
decomposition of scalar functions: Yi m (u)) = Yi m (uS), through the relation 
d l m0 (8) = [(/ — m)\/{l + m)!] 1 / 2 P i m (cos^). Any spin n function n G(u) on the 
sphere is thus uniquely given as a linear combination of spin n spherical har- 
monics : n G(uj) = Z)«eNS|m|<z nGim nYi m {oj) (inverse transform), for the spin- 
weighted spherical harmonics coefficients n Gi m = J S 2 dfl n Y^(a>)(jr(a;) (direct 
transform), with I > \n\, and \m\ < I. 

Finally, spin nil functions may be defined from spin n functions through the 
action of the so-called spin raising and lowering operators [T5p6] . The action 
of the spin raising S and lowering S operators on a spin n function n G, giving 
spin n + 1 and n — 1 functions respectively, is defined as 
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[3 n G] (6,<p) 



- sin™ 6 



d id 
+ 



86 sin 6 dip 



sin" n 9 n G 



•<p) (5) 



and 



3»G 



^ (96* sin (9<y9 / 



sin n # „G 



(6) 



with, under rotation by xo i n the tangent plane at cj = (9,<p): [8 n G]'(w) 



-i(n+l)xo 



[3 n G](w) and [3 „G]'(w) = e 



-i(n- l)xo 



9 n G](cj). In these terms, the 



spin-weighted spherical harmonics of spin n are related to spin- weighted spher- 
ical harmonics of spins n + 1 and n — 1 through the following relations: 



[3 n F im ] (c) = [(Z - n) (Z + n + 1)] 1/2 n+l Y lm (a) 



(7) 



and 
also implying 



9 nYlm 



(u) = -[(l + n) (Z-n + 1)] 



1/2 



.iFz m (a;) 



99 „H 



(a;) = - (Z - n) (Z + n + 1) n Y lm (to) . 



(9) 



The corresponding direct relation between the spin-weighted spherical har- 
monics of spin n and scalar spherical harmonics reads: 



(J-n) 
(Z + n) 



1/2 



[S n ^m] (w) 



(10) 



for < n < /, and 



(Z + n)! 



.1/2 



(Z-n)! 



/m 



for — Z < n < 0. 



These relations between spin-weighted and scalar spherical harmonics are ex- 
plicitly used in § [3] for the development of a fast direct spin ±2 spherical 
harmonics transforms algorithm. 



3 Fast spin ±2 transforms algorithm 



In this section, we define and implement a fast and exact algorithm for the 
computation of the spin ±2 spherical harmonics transforms of band-limited 
functions on equi-angular pixelizations on the sphere. The algorithm is based 
on the relations between spin-weighted and scalar spherical harmonics estab- 
lished in the previous section. 
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3.1 Pixelizations and existing 0(L 3 ) algorithms 

A 2L x 2L equi-angular pixelization in spherical coordinates (9, <p) is defined on 
points Uij = (6i, cpj) for < i, j < 2L — 1, with a uniform discretization of the 
coordinates: A9 = 9 i+ i — 9i = tt/2L and A<p = <pj + x — (pj = 2n/2L. The specific 
choice 9 = vr/4L, and (po = is considered in the following implementations. 
It gives 9i = (2i + l)ir/4L and (pj = 2jir/2L, and excludes the poles of the 
sampling, which can be convenient for numerical reasons. The pixels centers 
are identified with the sampling points u>ij defined here above. The pixels edges 
are identified by meridians shifted by A9/2 = n/AL, and parallels shifted by 
A(p/2 = 27r/4L relative to cjy. The poles therefore appear as pixels corners. In 
the next paragraphs, we analyze the properties of equi-angular pixelizations 
which are of interest in the implementation of scalar and spin ±2 spherical 
harmonics transforms. These properties are discussed in comparison with the 
HEALPix pixelization^ (Hierarchical Equal Area iso-Latitude Pixelization) 
[T8] . and the GLESP pixelizatiorO] (Gauss-Legendre Sky Pixelization) [I~9||2"0] . 
which are widely used in astrophysics and cosmology. 

Firstly, we discuss the asymptotic complexity for the computation of scalar 
and spin ±2 spherical harmonics transforms. Let us consider band-limited 
functions n G(uj) on the sphere with band limit L, defined through the fol- 
lowing condition on their scalar (n = 0) or spin-weighted (n 7^ 0) spherical 
harmonics coefficients: n Gi m = for / > L. For a signal with band-limit L, 
the a priori complexity associated with the naive computation of the direct 
scalar spherical harmonics transform integral on the sphere through simple 
discretization, i.e. a quadrature, for all (l,m) with \m\ < I < L, is naturally 
of order 0(L 4 ). And the a priori complexity associated with the naive com- 
putation of the direct spin ±2 spherical harmonics transforms integrals on the 
sphere through simple quadrature, for all (l,m) with I > 2, and \m\ < I < L, 
is also naturally of order (9(L 4 ). The same complexity naturally applies to the 
corresponding inverse scalar or spin ±2 transforms. We consider fine samplings 
corresponding to megapixels maps on the sphere. In particular, the WMAP 
experiment currently provides all-sky maps of around three megapixels. For 
such a fine sampling defining a band limit around L ~ 10 3 , the typical compu- 
tation time for (2L) 2 multiplications and (2L) 2 additions of double-precision 
numbers is of order of 0.03 seconds on a standard 2.2 GHz Intel Pentium 
Xeon CPU. We take this value as a fair estimation of the computation time 
required for one integration for given (l, m), or one summation for given (9, if), 
with an associated 0(L 2 ) asymptotic complexity. Consequently, scalar or spin 
±2 spherical harmonics transforms, with an asymptotic complexity of order 
0(L A ), typically take several days at that band limit L ~ 10 3 on a single 



1 http://healpix.jpl.nasa.gov/ 

2 http : / / www .glesp.nbi.dk/ 
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standard computer. Considering the analysis of a large number of signals or 
simulations may become difficultly affordable in terms of computation times, 
a fortiori in the perspective of forthcoming experiments with improved reso- 
lution on the sky, such as the Planck satellite experiment which will release 
all-sky maps of around fifty megapixels. 

The development of a fast and exact algorithm is therefore of great interest for 
the CMB analysis. The technique of separation of variables in the scalar or spin 
±2 spherical harmonics into the associated Legendre polynomials P/™(cos#) 
or the Wigner (i-functions d l m2 (6), and the complex exponentials e imip allows 
to decompose the transform as successive transforms in ip and 9 [2TH22] . It 
naturally reduces the asymptotic complexity for the direct and inverse scalar 
and spin ±2 spherical harmonics transforms to 0(L 3 ). It can be implemented 
on any iso-latitude pixelization. Many pixelization schemes have been con- 
sidered on the sphere which satisfy this requirement. It is the case for the 
equi-angular, HEALPix, and GLESP pixelizations. The algorithms existing 
on HEALPix or GLESP pixelizations are indeed based on this technique. As 
discussed in the next subsection, the asymptotic complexity may be further 
reduced on equi-angular pixelizations. 

Secondly, we discuss the precision of the computation. A sampling theorem 
exists on equi-angular pixelizations on the sphere, which represents a gener- 
alization of the Nyquist-Shannon theorem on the line. The sampling theorem 
states that the scalar spherical harmonics coefficients of a band-limited func- 
tion on the sphere may be computed exactly up to a band limit L, through 
a 2L x 2L equi-angular sampling, as a finite weighted sum, i.e. a quadra- 
ture, of the sampled values of that function (23]. The weights are defined from 
the structure of the Legendre polynomials P/(cos0) on the interval [0,7r]. A 
Gaussian quadrature rule for the exact computation of spherical harmonics 
coefficients of band-limited functions also exists on GLESP pixelizations. The 
HEALPix implementation of the scalar and spin ±2 spherical harmonics trans- 
forms achieves a very good precision thanks to an iteration process, but it is 
only approximate from the theoretical point of view as no sampling theorem 
is established on such pixelizations. 

Thirdly, we comment on the notion of pixel window function. On equi-angular 
pixelizations, the area A of pixels varies considerably with the co-latitude, 
from small pixels close to the poles, to larger pixels around the equator: 
A(uij) ~ sin OiAOAip. This is a major difference with the HEALPix pixeliza- 
tion which defines equal-area pixels, or the GLESP pixelization which defines 
nearly equal-area pixels. The constant area of pixels is an important property 
allowing the definition of a pixel window function associated to a given pix- 
elization at a given resolution. The main interest of this concept is to apply a 
low-pass filtering to the signal, implementing the fact that the pixelized signal 
is smoothed by integration over the pixel area. The corresponding window 
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function depends on the pixelization structure and resolution. The procedure 
of pixelization is approximated to a correlation of the signal with an axisym- 
metric beam, and therefore strongly relies on the assumption of equal-area 
pixels. We do not consider here the generalization of this concept on equi- 
angular pixelizations, where the pixel area varies drastically over the surface 
of the sphere. We only consider signals with band limit L on a 2L x 2L equi- 
angular sampling. In this case, for an application such as the downsampling, 
the spherical harmonics coefficients of a signal can be computed exactly thanks 
to the sampling theorem, and truncated at the desired band limit. In that re- 
spect at least, the use of the pixel window function can be avoided. 

Let us finally emphasize that each pixelization scheme (equi-angular, HEALPix, 
GLESP, ...) may provide specific advantages. The new algorithm proposed in 
the next subsection on equi-angular pixelizations is exact and has an asymp- 
totic complexity of order 0(L 2 log 2 L). But pixelizations with equal-area pixels 
represent an advantage when dealing with noisy data |18j . Our algorithm is 
therefore to be understood as a simple alternative to the existing algorithms. 
A more detailed comparison of the various algorithms is out of the scope of 
the present work. 



3.2 New exact 0(L 2 \og 2 2 L) algorithm 

We recall the following derivative relation on the associated Legendre polyno- 
mials [13], 



under the convention that P™ is defined to be zero for I < \m\. Through this 
relation, the derivative relations ffTU]) and (jlip between the spin ±2 spherical 
harmonics ±2^m an d the scalar spherical harmonics Yi m may be turned into 
a simple expression of ±2^m as linear combinations without derivatives of 
Yimi y«-i)mj an d y«_2)m- Notice that the same recurrence procedure is used 
in a different context in [21], in order to express spin n spherical harmonics 
nXirn, for any n with < \n\ < I, as linear combinations of scalar spherical 
harmonics. Through the recurrence relation on I satisfied by the associated 
Legendre polynomials of given m, 

(I -m)PT (cos 6) = (21- 1) cos (cos 0)-(Z + m- l)P£ 2 (cos0), (13) 

the Yn2)m term in the quoted linear combination for ±2^m may be cancelled. 
We finally obtain the following expression of ±2^m as a linear combination of 




(12) 
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the scalar spherical harmonics Y\ m and Yi 



(l-l)m- 



±2 



r lm (e,<p) 



(l-2)\ 



1/2 



.(i + 2)!. 

for Z > 2 and \m\ < I, and with the functional coefficients 



(14) 



4m) W 

(f ) 



2m 2 -Z(Z + l) , cot 9 ,„ , - 

; - =f 2m (Z — 1) — — + Z (Z - 1) cot 2 9 



sin 2 6 



sin 6 1 



2/ 



1/2 



m 



+ 



cot 9 s 



sin 9 sin 9 



(15) 



This relation holds once more under the convention that Y\ m is defined to be 
zero for Z < \m\. 



Consequently, the direct spin-weighted spherical harmonics transform of a 
spin ±2 function ± 2 G may be written as a linear combination of direct scalar 
spherical harmonics transforms for three associated scalar functions. Indeed, if 
the associated functions are defined by G^(9,ip) = (cot p 9/ sin 9 9)± 2 G(9, (p), 
for p, q G N and p + q = 2, one gets from relation (fT4l) : 



±2Glm 



[(1-2)!" 


1/2 (- 


L(l + 2)!_ 


i 2 



2Z + 1 
2Z~ 



11/2 



m 



+Z (Z - 1) G(^ m =f 2m (Z — 1) G« /m + [2m 2 - Z (Z + 1)] G^) lm } , 

(16) 



with Z > 2 and |m| < Z. The relation fj!4H also implies that the inverse spin- 
weighted transform of a set of spin ±2 coefficients ±iG\ m (with ±267™ = for 
Z > L) may be written as a sum of three inverse scalar spherical harmonics 
transforms: 

G (9, <p) = -±-A (9, <p) + ^B (9, <p) + cot 2 9C (9, <p) , (17) 
sin 9 sm 9 

with the scalar functions A, B, and C identified as follows by their scalar 
spherical harmonics coefficients: 
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A 



lm 



" (*-2)i 
(l + 2)\ 



1/2 



2m 2 -1(1 + 1) 



±2^m 



5, 



±2m 
'{1-2 



(/-!)! (2/ + 3) 



.1/2 



lm 



(l + 3)\ (2/+ 1 

1/2 

[ T 2m (/ - 1)] ±2 G lm 



2 2 

— m 



±2G ! (z + i) rn 



+2 



(J + 2)! 

(1-1)1 (21 + 3) 



(Z + 3)!(2Z + 1) 



m 



1/2 



(71 



lm 



(/-2)! 
(Z + 2)! 



1/2 



[I (l-l)}± 2 G lm . 



As discussed, for functions band-limited at L, the separation of variables al- 
lows to compute the direct and inverse scalar spherical harmonics transforms 
in 0(L 3 ) operations. However, a faster algorithm was developed by Driscoll 
and Healy on 2L x 2L equi-angular pixelizations on the sphere for the scalar 
spherical harmonics transforms [23]. The Fourier transforms in e irwf are com- 
puted in 0(L\og 2 L) operations for each 9 through standard Cooley-Tukey 
fast Fourier transforms. The algorithm also explicitly takes advantage of the 
recurrence relation in I on the associated Legendre polynomials P/"(cos#) to 
compute the direct associated Legendre transforms in 0(L\og\L) operations 
for each m. In these terms, the direct and inverse scalar spherical harmon- 
ics transforms are computed in 0(L 2 log 2 L) operations. The computation is 
theoretically exact thanks to the sampling theorem on equi-angular pixeliza- 
tions. Corresponding stable numerical implementations exist in the Sphar- 
monicKit package [25ll26f^1 . Through the relations (fT6l) and (IT7j) . the spin- 
weighted spherical harmonics transform of a band-limited spin ±2 function 
with band limit L may consequently also be computed exactly on a 21 x 2L 
equi-angular pixelization on the sphere from the Driscoll and Healy fast scalar 
spherical harmonics transforms, and with the same asymptotic complexity of 
order 0(L 2 log 2 L). In terms of our previous intuitive estimations, we recall 
that an 0(L 2 ) scalar product requires the order of 0.03 seconds on a standard 
2.2 GHz Intel Pentium Xeon CPU, at band limits around L ~ 10 3 . When com- 
pared to the a priori 0(L A ) asymptotic complexity, the 0(L 2 log 2 , L) scalar 
and spin ±2 spherical harmonics transforms algorithms consequently reduce 
computation times from days to seconds for the fine samplings considered. 

Let us remark that a recurrence relation was proposed in [27] in order to com- 
pute spin n spherical harmonics transforms from scalar spherical harmonics 
transforms. However, the proposed relation explicitly relates „Yj m with nl :iYi m , 
„ T iY(j-i) m , and „ T iY(i + i) m . The term „ T iF( i+ i) m increases the band limit of the 
functions to be analyzed to L + 2 after the 2-steps recurrence leading from 



http:/ /www. cs.dartmouth.edu/~geelong/sphere/ 
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spin ±2 to scalar spherical harmonics. On 2L x 2L equi-angular pixelizations, 
the SpharmonicKit package is technically limited to consider coefficients lower 
than L, and numerical errors will occur due to the absence of consideration of 
the coefficients at I = L and I = L + 1. No such issue occurs from the relation 
(fl4|) here above, which preserves the band limit L for the associated scalar 
functions. 



3.3 Numerical implementation 

We here report the computation times and memory requirements for the nu- 
merical implementation of the algorithm at band limits up to L = 1024, and 
briefly discuss the issue of the numerical stability of the implementation. The 
implementation is directly based on the fast scalar spherical harmonics trans- 
form proposed by Driscoll and Healy and implemented in the SpharmonicKit 
package. Computations are performed on a 2.20 GHz Intel Pentium Xeon CPU 
with 2 Gb of RAM memory. Random band-limited test-functions are consid- 
ered. Without loss of generality, these test-functions are defined through their 
spin-weighted spherical harmonics coefficients ±2^ m , with \m\ < I < L, and 
I > 2, with independent real and imaginary parts uniformly distributed in the 
interval [— 1,+1]. The inverse and direct spin-weighted spherical harmonics 
transforms are successively computed, giving numerical coefficients n Hi m . 

The computation times given in Table [1] for the direct and inverse spin ±2 
transforms are averages over 5 random test-functions. They range between 
1.0 x 10" 1 seconds for L = 128 and 2.2 x 10 1 seconds for L = 1024. The 
equality of computation times for the positive and negative spins is an evident 
consequence of the similarity of the ±2 cases in the relation ffl4|) . The case 
n = corresponds to the scalar spherical harmonics transform, and is added 
for comparison. The related values range between 2.7 x 10~ 2 seconds for L = 
128 and 6.5 seconds for L = 1024. To summarize, computation times are 
of the order of seconds for a band limit L = 1024, in agreement with our 
previous intuitive estimations. Both for the direct and inverse transforms, the 
evolution of the values reported as a function of the band limit also supports 
the 0(L 2 \ogl L) behavior of the related asymptotic complexity, as illustrated 
in figured] in comparison with an 0(L 3 ) slope. The ratio of computation times 
for the cases n = ±2 and n = also reflects the simple fact that three scalar 
transforms are computed for each spin ±2 transform. 

In the present implementation based on the SpharmonicKit package, the re- 
quired associated Legendre polynomials P™(cos0) are pre-calculated once for 
all values of I, 8, and m, and stored in RAM memory. The pre-computation 
time itself is of order 0(L 3 ) through the use of a recurrence relation in I on 
the associated Legendre polynomials. This pre-computation is by definition 
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Table 1 

Computation times for n = and n = ±2 spherical harmonics transforms measured 
on a 2.20 GHz Intel Pentium Xeon CPU with 2 Gb of RAM memory. Times asso- 
ciated with the direct transforms are listed above the corresponding times for the 
inverse transforms. 
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Figure 1. Evolution of computation times t displayed as log 2 1 — log 2 L for the direct 
(left) and inverse (right) spin-weighted spherical harmonics transforms of spins n = 
(continuous line), n = +2 (continuous line), and n = —2 (dashed line). Computation 
times are measured in seconds on a 2.20 GHz Intel Pentium Xeon CPU with 2 
Gb of RAM memory, and reported for the band limits L G {128,256,512,1024}. 
The 0(L 2 log 2 L) asymptotic complexity is clearly illustrated when compared to an 
£>(L 3 ) slope (dotted line). 

not taken into account in the reported computation times, which consequently 
remain of order 0(L 2 log 2 . L). The number of real values of associated Legen- 
dre polynomials P; m (cos^) stored in RAM memory for all I, 8, and m is also 
of order 0(L 3 ). The overall memory requirements allowing the direct and in- 
verse transforms with the present numerical implementation correspondingly 
increase from 5.6 Mb for L = 128, to 32 Mb for L = 256, 220 Mb for L = 512, 
and 1.2 Gb for L = 1024. These memory requirements are easily accessible on 
a single standard computer. 
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Table 2 

Errors are measured on a 2.20 GHz Intel Pentium Xeon CPU with 2 Gb of RAM 
memory. Absolute errors after inverse and direct transforms are listed above the 
corresponding relative errors. 



The absolute and relative numerical errors are defined as max; m I nGi m — 
n H lm \ and max l])7l |( n G lm - n Hi m )/ n Gi m \ respectively, where | • | here denotes 
the complex norm, and n G {0, ±2}. The numerical errors associated with the 
0(L 2 logg L) spin- weighted spherical harmonics transforms given in Table [2] 
are averages for transforms over 5 random band-limited test-functions. Abso- 
lute and relative errors do not exceed the order of 8.4 x 10~ 9 and 4.2 x 10~' 
respectively for band limits up to L = 1024. The Oil? log^ L) implementation 
of the spin ±2 spherical harmonics transforms is therefore stable for band lim- 
its up to L = 1024. The numerical stability of the algorithm might also have 
been inferred from the corresponding stability of the Driscoll and Healy fast 
direct scalar spherical harmonics transform algorithm, tested for band limits 
up to L = 1024 [25 26J. The only potential source of instability is related to the 
multiplication factor cot p 9/ sin 9 9 defining the scalar functions associated with 
a spin ±2 function in the computation of a spin-weighted spherical harmonics 
transform from the relation (|T4l) . Each such factor indeed corresponds to a 
division by sin 2 9, which induces multiplications by numbers of the order of L 2 
around the poles 9 = {0,7r}, where L is the band limit considered. However 
such operations could only produce numerical instabilities for very high band 
limits, and obviously remain completely safe at L = 1024. 



4 Application in cosmology 

In this section, we illustrate the interest of the algorithm presented in the 
previous section in the context of the analysis of CMB polarization data. 
The discussion is based on the following introductory papers [5 6 7 8 .9 TfT0] and 
reviews [28,29,30 31] relative to the CMB polarization analysis. 
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4-1 Stokes parameters 



The CMB is observed in each direction uo = (9, p) of the sky as an incoming 
radial radiation, to which is associated a transverse electromagnetic field thus 
lying in the tangent plane to the sphere at the point considered. In that plane, 
we consider the basis (eg, e v ) with eg pointing in the direction of increasing 9 
along each meridian, and e v in the direction of increasing p along each paral- 
lel. In this so-called linear polarization basis, nearly monochromatic radiation 
around a frequency u r may be decomposed as an electric field with compo- 
nents E e (u r , t) = TZe[ee(t)e~ iuJrt ] and E v (u r , t) = fte[^(t)e" w ]. The complex 
amplitudes Sg(t) and s v (t) slowly vary in time relatively to the timescale set by 
the wave period. The intensity matrix I associated with the radiation simply 
reads as the time average of the electric field rank 2 tensor [e*(t)ej(t)]ii <8> ij, 
for i,j e {9, ip} [5][T0] . It thus naturally decomposes on the 2x2 matrix ba- 
sis formed by the identity matrix cr = I, and the well-known Pauli matrices 
cr 2 , 03), as I = [Ia + Uo\ + Va 2 + Qa 3 }/2. The constants J, U, V, and 
Q define the four real Stokes parameters [32J associated with the radiation: 

1 = (|e*(*)| 2 + M*)| 2 ) > Q = (\ £ e( t )\ 2 -K( t )\ 2 )^ u = ( E *e( t )^^( t )+^e(t)e* v (t)), 
V = i(eQ(t)e tp (t) — Eg(t)e* (t)). The brackets (■) denote time averaging. If the 
two components eg(t) and e v (t) are correlated, the radiation is said to be po- 
larized. The positive parameter / may be identified with the overall intensity 
of radiation, while Q and U identify with the linear polarizations, and V with 
the circular polarization. Unpolarized radiation, or natural light, is therefore 
characterized by Q = U = V = 0. 

As functions on the sphere, Q(u>), U(uS) and V(uj) have different be- 

haviors both under parity, i.e. global inversion (•") of the coordinates, and 
under local rotations (•') of the basis vectors (e^e^) in the tangent plane at 
u — (8, f). A global inversion of the right-handed three-dimensional Cartesian 
coordinate system (o, ox, oy, oz) centered on the unit sphere induces the fol- 
lowing modification of Cartesian coordinates: (x", y" , z") = (—x, —y, —z). The 
spherical coordinates uj = (6, ip) of a given point on S 2 change according to 
oj" = (8", p") = (n—O, n+tp). Locally in the tangent plane, the global inversion 
also implies an inversion of the basis vector eg: (e'g, e") = (—eg, e^). The Stokes 
parameters I and Q have even parity, I"(uj") = I(uj) and Q"(oj") = Q(ui), 
while U and V have odd parity, U"(u") = -U(u) and V"(u") = -V(u). Un- 
der local rotations of the basis vectors (eg, e v ) by an angle Xo> the coordinates 
e = (sgjEip) of vectors in the tangent plane transform through e' = r xo ■ e, 
for the standard rotation matrix r xo , with entries = = cos^o an d 
= — r^h = sinxo- The Stokes parameters / and V are invariant while Q 
and U are mixed by local rotations. Equivalently, one may also rewrite the 
intensity matrix as 

I = i [Jo-o + Va 2 + (Q + iU) a + + (Q- 1U) <r_] , (19) 
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with cr± = ((T3 =Fio"i)/2. The Pauli matrices transform as a'^ = r Xo -cr^-r^, for 
/x = {0,l,2,3}. The matrices o"o and o"2 are thus invariant, while a± transform 
as a± = e T2lXa a±. Consequently the four Stokes parameters are associated with 
spin functions on the sphere. The intensity I{ui) and the circular polarization 
parameter V(u) are scalar functions. The combinations (Q ± iU)(u) are spin 
±2 functions: (Q±iU)'(uj) = e T2%xo (Q±iU)(oj). Notice that under parity these 
two combinations transform in one another: (Q ± iU)"(u") = (Q =F iU)(ui) 



4-2 Angular power spectra 



It is assumed that the physics of the CMB is invariant under parity and under 
local rotations. It is therefore suitable to relate the observables I, Q, U, and V 
to invariant physical quantities. The intensity I(co) defines the CMB tempera- 
ture anisotropics T(u) and is indeed itself invariant under the transformations 
considered. As no circular polarization may arise from Thomson scattering, 
the CMB polarization is completely described in terms of the two linear polar- 
ization Stokes parameters Q and U. It is equivalently defined by their spin ±2 
combinations Q ± ill. Associated polarization components, real scalar func- 
tions on the sphere and parity eigenmodes, are naturally defined from Q ±iU 
in terms of the raising 9 and lowering 8 operators respectively given in (JSj) 
and ©. These components E(u) = -[8 2 (Q + iU)(u) + S 2 (Q - iU)(u)}/2, 
and B(uj) = i[Q 2 (Q + iU)(u) — <S 2 (Q — iU)(u>)]/2, have even and odd parities 
and are therefore referred to as electric and magnetic components respec- 
tively Let us consider the decomposition of the spin ±2 functions Q ±iU 
in spin-weighted spherical harmonics and the relations 2^zm = N^ 2 )Q 2 Yi m 
and _ 2 Y lm = N [l2) Q 2 Y lm , with N {12) = [(/ - 2)!/(Z + 2)!]^ irid uced from 
(TlO|) and (TTTj) . The application of the raising and lowering operators on this 

decomposition through the relations (JTj) to Q) gives E[ m = Ei m /N^ and 
Bi m = B lm /N(i2), where 

Elm = ~\ ( + 2{QTTu) lrn + _ 2 (Q^XJ) lm ) (20) 

and 

Bim = \ ( + 2{Q + iU) lm - _ 2 {Q^Tu) lm ) , (21) 

define the properly normalized real E(u) and B(u) components. These coef- 
ficients are explicitly invariant under local rotations. 

The random process from which the CMB radiation arises is assumed to be 
Gaussian and stationary. It is therefore completely characterized in terms of its 
temperature and polarization two-point correlation functions. The correspond- 
ing invariant angular power spectra are naturally those associated with the 
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temperature TT, the polarizations EE and BB, and the cross-correlation be- 
tween the temperature and electric polarization component TE: (Tn m ,Ti m ) = 

mm' ; 

Cf E 5u'5 mm '. These physical quantities are indeed invariant under local rota- 
tions and parity. The TB and EB cross-correlations are specifically excluded 
from the requirement of invariance under parity 

Notice that the E and B components of polarization not only define invariant 
physical angular power spectra, but they are also associated with different 
mechanisms of production of the radiation, corresponding to different theoret- 
ical cosmological models. Scalar primordial energy density perturbations only 
produce the E polarization component, while vector and tensor (i.e. gravity 
waves) perturbations produce both E and B polarization components. 

4-3 Numerical illustration 

Scalar and spin ±2 direct spherical harmonics transforms are required for 
the estimation of the CMB angular power spectra from the observables T, 
Q, and U [6]. The simulation of temperature and polarization maps from 
given theoretical angular power spectra requires the corresponding inverse 
transforms. We apply our algorithm to simulate CMB maps and angular power 
spectra, for illustration of its precision and speed performances. 

We start from the temperature and polarization spectra Cf T , C\ ■ , and Cf E 
defined by the concordance cosmological model which best fits the three-year 
WMAP data (the BB polarization spectrum is identically null: Cf B = ). 
These spectra are represented in figure[2]up to a band limit L = 1024. Spherical 
harmonics coefficients TJ m and E[ m are built up as the two marginal complex 
Gaussian realizations arising from a jointly Gaussian statistical distribution 
with variances Cf T and C EE , and a covariance Cj E . The T , Q, and U maps are 
then produced by inverse scalar and spin ±2 transforms, through the relations 
(EOD and (J2I1), with B lm = 0. 

From the maps obtained, we recompute spherical harmonics coefficients T( m , 
E' lm , and B' lm by direct scalar and spin ±2 transforms. Within the numeri- 
cal accuracy of the computer, the B polarization coefficients are identically 
null, in perfect agreement with the original data: B[ m = 0. We finally es- 
timate the temperature and polarization angular power spectra from those 
coefficients as: C? T > = Y? m =-i l*Ll7(21 + 1), Cf B ' = E l m =-i IKJ7(2/ + 1), 
CJ EI = EL-,nA/(2/ + l), and C? B > = E l m =-i \B' lm \ 2 /(2l+l) = 0. These 
estimators follow chi-square distributions with 21 + 1 degrees of freedom. For 
X G {TT, EE, BB,TE}, this induces a fractional uncertainty o c x,jCf = 

[2/(2/ + l)] 1 / 2 in the estimation. This cosmic variance is large at low I and 
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Figure 2. CMB temperature and polarization angular power spectra Cf T (top), 
C EE (bottom left), and Cf E (bottom right) of the cosmic microwave background up 
to a band limit L = 1024 and in (iK 2 . Inverse and direct transforms are successively 
performed through the exact 0{L 2 \og^L) scalar and spin ±2 spherical harmonics 
transforms on 2L x 2L equi-angular pixelizations on the sphere, in order to produce 
the estimated spectra (scattered points) from the original spectra of the concordance 
cosmological model (continuous lines). The original and estimated spectra coincide 
within the 3a uncertainty defined by the cosmic variance (grey region). 

small at high /. The figure [2] represents the good coincidence between the orig- 
inal and estimated spectra up to the corresponding uncertainty at each I. The 
computation time associated with the overall procedure is 150 seconds on a 
2.20 GHz Intel Pentium Xeon CPU with 2 Gb of RAM memory. In summary, 
this application illustrates the good precision and speed performances of our 
fast and exact algorithm, coherently with the results of tables [1] and [2j 



5 Conclusion 

In conclusion, we developed a fast, exact, and stable algorithm for the spin ±2 
spherical harmonics transforms of band-limited functions with band limit L on 
2L x 2L equi-angular pixelizations on the sphere. The algorithm is based the 
Driscoll and Healy fast scalar spherical harmonics transform algorithm. The 
exactness of the computation on equi-angular pixelizations relies on a sampling 
theorem. The associated asymptotic complexity is of order 0(L 2 log*. L). The 
algorithm is presented as an alternative to existing algorithms with an asymp- 
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totic complexity of order 0(L 3 ) on the HEALPix and GLESP pixelizations, 
which are widely used in the context of astrophysics and cosmology. 

The numerical implementation produced confirms the characteristics of the 
algorithm. Typical computation times for L = 1024 are of the order of seconds. 
We also illustrated the interest of the algorithm in the context of the analysis 
of CMB polarization data. 
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